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Satellite remote sensing of ocean color in dynamic coastal, inland, and nearshore waters is impeded by high var- 
iability in optical constituents, demands specialized atmospheric correction, and is limited by instrument sensi- 
tivity. To accurately detect dispersion of bio-optical properties, remote sensors require ample signal-to-noise 
ratio (SNR) to sense small variations in ocean color without saturating over bright pixels, an atmospheric correc- 
tion that can accommodate significant water-leaving radiance in the near infrared (NIR), and spatial and tempo- 
ral resolution that coincides with the scales of variability in the environment. Several current and historic 
space-borne sensors have met these requirements with success in the open ocean, but are not optimized for 
highly red-reflective and heterogeneous waters such as those found near river outflows or in the presence of sed- 
iment resuspension. Here we apply analytical approaches for determining optimal spatial resolution, dominant 
spatial scales of variability (‘'patches"), and proportions of patch variability that can be resolved from four river 
plumes around the world between 2008 and 2011. An offshore region in the Sargasso Sea is analyzed for com- 
parison. A method is presented for processing Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua 
and Terra imagery including cloud detection, stray light masking, faulty detector avoidance, and dynamic aerosol 
correction using short-wave- and near-infrared wavebands in extremely turbid regions which pose distinct op- 
tical and technical challenges. Results show that a pixel size of -520 m or smaller is generally required to resolve 
spatial heterogeneity in ocean color and total suspended materials in river plumes. Optimal pixel size increases 
with distance from shore to -630 m in nearshore regions, -750 m on the continental shelf, and -1350 m in the 
open ocean. Greater than 90% of the optical variability within plume regions is resolvable with 500 m resolution, 
and small, but significant, differences were found between peak and nadir river flow periods in terms of optimal 
resolution and resolvable proportion of variability. 

© 2013 Elsevier Inc. All rights reserved. 


1. Introduction and theory 

In measuring the length of a particular stretch of convoluted 
coastline with many small bays and inlets, very different results are 
reached when measuring by the centimeter, meter, or kilometer. If 
the exercise were repeated at a comparatively straight and feature- 
less coastline, these results may converge. Similarly, the composition, 
distribution, and diffusion of the components of seawater that give it 
color— chromophoric dissolved material, suspended minerals, phyto- 
plankton, etc.— will exhibit patchiness on scales from microns to 
many kilometers. Any empirical approach to the question of finding 
the “best” resolution for ocean viewing is limited by the resolution 
at which observations are made. In this study, we address the issue 
using what is currently the most readily available, highest resolution 
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ocean color imagery with global coverage: 250 m MODIS data for the 
645 nm band. Using thousands of scenes from environments as 
dynamic as turbid river plumes after tropical storm flooding events 
to the oligotrophic Sargasso Sea, statistical and analytic tools are 
adapted to quantify interpixel variability from instrument noise, 
atmospheric correction, and ocean color itself. The spatial scales of 
variability (“patchiness”) are identified in the process, as well as the 
upper limits of the pixel size required to capture dispersion of 
suspended sediments driving ocean color, and the degree of variabil- 
ity that can be resolved given these limitations. 

1.1. Background 

A recent study by Lee, Hu, Arnone, and Liu (2012) demonstrated 
that satellite ocean color products derived from low resolution imag- 
ery are not necessarily the arithmetic or geometric average of those 
derived from high resolution imagery over the same scene. As a 
result, low resolution imageiy in environments with steep gradients 
such as eddies, fronts, phytoplankton blooms, and river plumes can 
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underestimate the concentration of retrieved optical properties and 
biogeochemical constituents. Using comparisons of the absorption 
coefficient retrieved from high resolution (300 m) and low resolution 
(1.2 km) Medium Resolution Imaging Spectrometer (MERIS) imagery 
over a large (-21,800 km 2 ) eddy in the Southern Ocean, they found 
that absorption was underestimated in the low resolution product 
by ~5-15%— particularly along the strongest spatial gradients in the 
absorption coefficient. By extension, they conclude that global esti- 
mates of satellite-retrieved products (i.e. chlorophyll concentration, 
Chi) have historically been underestimated due to the relatively 
coarse resolution of the current suite of satellite sensors. 

Patchiness in optical constituents of the ocean is present at all spa- 
tial scales. Physical and biological oceanographic processes also occur 
at spatial scales from sub-millimeter to thousands of kilometers 
(e.g. Fig. 6 in (Dickey, Lewis, & Chang, 2006)). Regions with strong 
gradients in Chi— such as a bloom— or total suspended material 
(TSM)— such as a river plume— will be inherently patchier than 
oligotrophic, low mineral content waters far from shore. The design 
of remote sensing satellites must be sensitive to the spatial scales of 
constituent variability in order to accurately assess their concentra- 
tion and distribution. As an example, in planning for the GEOstation- 
ary Coastal and Air Pollution Events satellite (GEO-CAPE) (Fishman et 
al., 2012), the technological, financial, and logistical costs of building a 
high spatial resolution sensor must be balanced against its capability 
to resolve natural variability and uncertainties in products retrieved 
in patchy environments. Nevertheless, high resolution (e.g. <1 km) 
may not be necessary for large areas of open ocean. Addressing the 
ocean science questions for GEO-CAPE (http://geo-cape.larc.nasa. 
gov/ocean.html) is a particular challenge with respect to spatial 
resolution because its high geostationary orbit (35,786 km altitude 
over -95° W) will require a more massive telescope and longer 
“stare interval” for any given target, which in turn limits coverage. Es- 
timating spatial resolution requirements for various environments 
from dynamic river plumes to the open ocean— as we endeavor to do 
here— can facilitate optimizing the revisit time over each target area, 
thereby increasing temporal resolution as well as spatial coverage. 

1.2. Tools for assessing GSD 

The optimal ground space distance (GSD) is the spatial resolution 
below which spatial heterogeneity in ocean color— or the underlying 
bio-optical constituents such as Chi or TSM— cannot be resolved. The- 
oretically, this decreases from quiescent regions offshore to coastal 
environments and turbid river plume regions. The latter regions are 
characterized by strong seasonal and episodic influx of TSM (mainly 
non-algal sediments), and present the additional challenge of being 
environments in which standard approaches to satellite imagery pro- 
cessing are often inadequate (Ruddick, Ovidio, & Rijkeboer, 2000). 

Ideally, any interpixel variability of adjacent pixels in ocean color 
products represents real differences in constituent concentrations. In 
such a case, finding the GSD for a given region would be a straightforward 
matter of reducing spatial resolution in an image (i.e. increasing GSD by 
averaging adjacent pixels) until variability in adjacent pixels is no longer 
negligible, beyond which the capacity for gauging dispersion in the con- 
stituent concentration is reduced. For example, if variability of adjacent 
pixels is significant among 250 m pixels, then larger pixels (e.g. 500 m) 
would overlook these differences. In fact, some of the interpixel variabil- 
ity in satellite images derives from non-geophysical influences such as 
the signal to noise ratio (SNR) of the satellite detector, artifacts arising 
from processing the image for atmospheric correction, and uncertainties 
in the algorithm used to retrieve the biogeochemical constituents from 
above water reflectances. Therefore, in order to determine optimal GSD 
in naturally noisy imagery, a threshold of noise above which differences 
in interpixel variability result from true differences in ocean color must 
be determined. To eliminate the complication of uncertainties introduced 
by ocean color algorithms, we generally limit our analysis here to the 


ocean color signal itself. Specifically, we will focus on the remote sensing 
reflectance at 645 nrn, R rs ('645 ). This waveband was chosen for three rea- 
sons, 1 ) it is available on both MODIS Aqua and Terra instruments and 
therefore has readily accessible global coverage together with up-to- 
date calibration and on-orbit signal-to-noise information, 2) it is mea- 
sured at relatively high resolution (250 m), and 3) it is an excellent 
proxy for TSM (e.g. (Kirk, 1994; Miller & McKee, 2004; Ondrusek et al., 
2012)) although regional optimization is generally required in MODIS 
band 1 (645 nm) algorithms for TSM as absorption features such as Chro- 
mophoric Dissolved Organic Material (CDOM) and phytoplankton pig- 
ments chlorophyll a and phycocyanin can significantly impact sea 
surface reflectance in this band. The 645 nm band in MODIS was original- 
ly designed as a land band, and therefore has a much lower SNR than the 
MODIS ocean bands (160.5 ± 7.2 at typical 645 nm oceanic radiance of 
1.72 ± 0.21 mWcm“ 2 pm -1 sr _1 versus -2000 in ocean bands (Flu et 
al., 2012 ) ). As we demonstrate in our results below, this introduces signif- 
icant noise in band 1 offshore, but in turbid, highly reflective river plumes 
such as those studied here, radiances are often five times greater, thereby 
boosting SNR to - 500 or higher. The 645 nm band on MODIS is also quite 
broad at -50 nm, but this should not impact interpixel variability for the 
purpose of GSD evaluation, nor is it likely to significantly impact its utility 
in deriving TSM since ocean color in this region of the spectrum is dom- 
inated by scattering— i.e. the single scattering albedo (the ratio of scatter- 
ing to total attenuation) approaches 1 in most environments (Babin, 
Morel, Fournier-Sicre, Fell, & Stramsld, 2003). 

Bissett et al. (2004) studied the issue of optimal GSD using aircraft 
imagery collected over the LEO-15 study site off New Jersey and 
found the random noise component of the surface reflectance signal 
by calculating the standard deviation between adjacent pixels in a re- 
gion within the image known to be relatively homogeneous in its 
water-borne optical properties, using this observed variability over 
a quiescent region as their threshold (cr t ). Assuming the noise to be 
linearly additive (constant, regardless of the magnitude of the signal), 
a transect of virtual stations was mapped onto the image in a nearby 
region in which optimal GSD was sought. At each station, an analysis 
was conducted in which the interpixel variability of pixels immedi- 
ately adjacent to the station is calculated (i.e. the standard deviation 
of a 3 x 3 array of pixels centered on the ith station, a,). If a, was 
higher than the threshold of random noise, a t , then it was determined 
that true differences existed in the optical properties of the water col- 
umn between these adjacent pixels, and a coarser GSD would miss 
this variability. If not, the array was increased to 5 x 5, 7 x 7, etc. 
until that condition was met. The resolution of the last aggregate 
array for which a, < a t determined the optimal resolution, or GSD. 
Based on a virtual transect superimposed on a single aircraft image 
collected by the PH1LLS 2 instrument at 9 m resolution at the 
LEO-15 study site, they found that near shore (<10 km from the 
coast) the optimal GSD was -100-200 m, and increased to -2000- 
6000 m between 10 and 50 km from shore. 

Here we pursue an approach similar to Bissett et al. (2004), but ex- 
pand it to include thousands of satellite images collected over several 
optically diverse regions and seasons. For practical reasons described 
in more detail below, the methods are modified in our analysis to ac- 
commodate the use of MODIS Aqua and Terra satellite imagery. Rather 
than using an optically quiet region of an image to estimate a t — as in 
(Bissett et al., 2004)— in this investigation, the noise threshold is deter- 
mined explicitly by combining on-orbit SNR for the MODIS band 1 with 
noise introduced by atmospheric correction. The optimal GSD in this 
study is established as the average between the size of the last pixel 
array for which interpixel variability is within the noise (a, < cj t ) and 
the size of the first pixel array for which variability exceeds the noise 
threshold (a ( > cr t ). A limitation in using MODIS QKM (nominally 
250 m) imagery is that it cannot extend to resolutions lower than the 
smallest aggregation of pixels— a 2x2 array— or an optimal GSD of 
375 m. In other words, recommendations made here regarding regions 
with the most stringent GSD requirements (i.e. plume regions) are 
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typically upper limits of pixel size, but both upper and lower bounds 
may apply nearshore, on the shelf, or offshore where ranges of optimal 
GSD are well above the 375 m limit. 

A statistical method often used to describe spatial variability in en- 
vironmental measurements is the semivariogram which builds on the 
intuitive premise that closer observations have a higher probability of 
being similar than distant observations. The shape of semivariance 
(y) with distance (or “lag” h, see Fig. 1) has been used to describe 
geographical patterns in the marine environment such as chlorophyll 
variability in the oceans (Doney, Glover, McCue, & Fuentes, 2003), to 
help establish spatial scales of variability in remotely sensed parame- 
ters, and to describe the proportion of variability captured at a given 
resolution (Atkinson & Curran, 1997; Davis, Kavanaugh, Letelier, 
Bissett, & Kohler, 2007). 

The semivariance function is found by calculating half the squared 
difference of all pixel pairs in an image at distance h from each 
other in any direction, over all possible distances. The example 
semivariogram in Fig. 1 shows idealized models for y (h) with 
descriptive characteristics. As distance h increases to the range, the 
semivariance first increases as more distant observations become 
less similar, and then plateaus at the sill— the distance at which differ- 
ences between observations are no longer a function of distance, but 
rather reflect the overall variability in the region (or sub-regional 
patch) being evaluated. The smallest distance in a satellite image is 
the same as the pixel size, so the nugget (i.e. the semivariance 
approaching the y-intercept) reflects the unresolved variability at a 
given resolution. The rate at which y rises to reach the sill (the 
shape of the curves in Fig. 1 ) can be viewed as an indication of inter- 
mediate scale patchiness, and nodes— deviations from the smooth, 
approximately hyperbolic slope— indicate spatial scales at which 
patchiness is a factor. Multiple sills may exist in a semivariogram 
which encompasses diversity at different spatial scales (e.g. dot- 
dashed line in Fig. 1). Typically, semivariograms assume an isotropic 
distribution of observations, so strong spatial gradients may distort 
results. If a gradient is known to exist a priori, it is also possible to 
calculate the non-isotropic semivariance; the semivariance along a 
particular direction believed to be free from significant, systematic 
gradients. Semivariograms are used in this study assuming isotropy 
within the smallest patches observed (typically 20-100 km)— and 
only within the plume regions— to establish spatial scales of patchi- 
ness, overall variability of a region, and the degree of variability 
captured by various resolutions of satellite imagery, while avoiding 
the strong gradient between the plume and open ocean waters 
which could distort results. 


-Sill (c,) 


Interpixel Distance (h) 

Fig. 1. Idealized semivariograms (dotted and dash-dotted lines) demonstrate spatial 
variability in nature tending to increase with distance between observations. The sill 
(ct) represents the total scale of variability of a region at a distance equal to the range 
(a) where variance plateaus. Increasing the size of the region to incorporate diversity 
not previously encompassed may cause variance to rise above one sill until a new sill is 
reached (dash-dotted line), demonstrating spatial scales of “patchiness”. If the 
semivariance is scaled to the sill in a given region, the partial sill (=1— nugget (Co)) is 
the proportion of variability that can be resolved as the distance approaches the y-inter- 
cept. This distance is equivalent to the minimum distance possible in a satellite image 
(i.e. the pixel size or resolution). 


c 

n 

Cu 


•' Sub-region Sill (Patch) 


ii; 

3 


Range (a) 


1.3. Uncertainly in surface reflectance 


The ability to determine optimal GSD from satellite imagery 
depends on establishing the interpixel variability in the reflectance 
signal above the variability caused by signal uncertainty or noise 
(i.e. Oj > CT t , where ct is the standard deviation of an array of pixels 
divided by the average: the coefficient of variability). For the spectral 
remote sensing reflectance (R rs (X); the water-leaving radiance scaled 
to the surface downwelling irradiance), uncertainty (A) can be 
expressed by propagation of error as: 


AR rs (\) = 


LwW 

EiW 


/ 


\AE d (\)J 


+ 


( AE d (\) \ 2 

\ EdW ) ’ 


( 1 ) 


where ^(A) is the spectral water-leaving radiance, E d (A ) is the spectral 
sea-surface irradiance, and \ is wavelength (spectral notation 
suppressed for brevity hereafter). Top of atmosphere radiance is the 
sum of radiance contributions from aerosols L a , glint L g , combined 
Rayleigh correction and Rayleigh-aerosol interactions L ra , white-caps 
L wc , and L w . Including the uncertainty in each of these elements as 
well as the atmospheric transmittance coefficient t yields: 

L t ± AL f = L a ± AL a + L g ± AL g + L m ± AL ra + L wc ± AL WC + t(L w ± AL W ), 

and again by propagation of error, 

jAL 2 t -AL 2 a -ALl-AL 2 m -ALl c 

AL t = V -L . (2) 


Combining Eqs. (1) and (2) yields total uncertainty in R rs : 


AR r , - R r 


AL 2 -AL 2 -AL 2 g -AL 2 a -AL 2 wc 
(f^-vv) 2 



(3) 


Given that the 645 nm MODIS band was designed for land use, its 
SNR ratio over water will be a major source of uncertainty, particular- 
ly in low turbidity pixels throughout our regions. Furthermore, as de- 
scribed in Section 2.1.5, atmospheric correction methods will depend 
on both near infra-red (NIR) and short-wave infrared (SW1R; 
1240 nm and 2130 nm for MODIS) bands. Due in large part to the 
low SNR in the SWIR bands (48.4 ± 3.9 and 30.8 ± 3.3, respectively 
(Hu et al., 2012)), atmospheric correction is likely to be another 
major source of uncertainty. Isolating the uncertainty in L t driven by 
instrument noise from the remaining atmospheric correction terms 
in Eq. (3) yields 

AR rs = AR“ + AR? S C , (4) 

where AR“ is uncertainty owing to the SNR in L t , and ARfif is uncer- 
tainty in the atmospheric correction. 

If uncertainty in R rs were solely due to instrument noise in the TOA 
radiance, where AL t = L t /SNR U , and SNR L[ is the signal to noise ratio 
in L f , Eq. (3) reduces to: 


AR " ~ Ks IlJnR,, 


(5) 


However, atmospheric correction (i.e. the remaining terms in Eq. (3)) 
introduces large uncertainty near shore and in turbid waters. Uncertain- 
ty in ocean color products deriving from atmospheric correction, includ- 
ing the SNR in NIR and SWIR bands, was estimated by the Pre-Aerosol, 
Clouds, and ocean Ecosystem (PACE) Science Definition Team in their 
report in October 2012 (http://decadal.gsfc.nasa.gov/pace.html). 
Interpolating between the 555 nm and 670 nm bands in that study 
(see their Fig. A-l), conservatively estimating SNR in NIR and SWIR 
bands from MODIS (-800 and -50, respectively (Hu et al., 2012)), and 
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converting to error in R rs from normalized water leaving radiance yields 
uncertainty in the 645 nm of approximately ARff =6.4x1 CP 5 sr^ 1 
when NIR bands are used for atmospheric correction, and AR?^ = 
8.0 x 10“ 4 sr* 1 when SWIR bands are used. 

From Eqs. (4) and (5), and a pixel-by-pixel assessment of which 
bands were used for atmospheric correction (Section 2.1.5), we esti- 
mate the non-biogeochemical variability in sea surface reflectance 
(cj t ) on a pixel-by-pixel basis for all images studied. 

2. Data, areas of interest, and methodology 

Prior studies of spatial resolution such as those mentioned above 
have focused on models and/or small sample sizes of imagery 
(Atkinson & Curran, 1997; Bissett et al., 2004; Davis et al., 2007; Lee et 
al., 2012; Ruddick et al„ 2012). Our goal is to adapt and generalize 
these approaches to cover plumes of major and minor rivers, nearby 
estuarine, coastal, and shelf waters, and regions far from land. 

2.3. Satellite imagery 

3076 MODIS level 1A (LI A) Aqua and Terra granules collected dur- 
ing peak and nadir flow periods over major river outflows (Amazon, 
Mississippi, Susquehanna, and Yangtze Rivers) between 2008 and 


2011 were downloaded from NASA's archive (http://oceancolor.gsfc. 
nasa.gov) and regionally extracted and processed using the SeaWiFS 
Data Analysis System (SeaDAS 6.4) to level 2 (L2) using a novel method- 
ology optimized for highly turbid zones, as described in Section 2.1.5 
and Fig. 3. In addition, 1104 LI A images from an optically quiescent re- 
gion in the Sargasso Sea bounded by 25°N and 35°N latitude and 
— 45°W and — 55°W longitude during May and November 2008-201 1 
were processed to L2 using standard operational methods for ocean 
color (Ahmad et al., 2010; Bailey, Franz, & Werdell, 2010; Gordon & 
Wang, 1994) L2 images were generated at 1000 m resolution (1KM), 
500 m resolution (HKM) and 250 m resolution (QKM), which is the na- 
tive, nadir, nominal resolution for the 645 nm band (band 1, MODIS) 
used here as an analog to TSM. Actual spatial resolution is a function 
of satellite viewing geometry, and is calculated for each image prior to 
GSD and semivariance analysis. QKM scenes with > 500 m native resolu- 
tion are rejected from analysis. 

Regions of interest are shown in greater detail in Fig. 2 together 
with example retrievals of R rs ( 645 /determined using the methods de- 
scribed below (Section 2.1.5). Sub-regions (plume, coastal/estuary/ 
nearshore, and continental shelf waters; black, red, and green boxes, 
respectively) are delineated in the figure, and the virtual transects 
and stations used in the analysis of optimal GSD in each sub-region 
are superimposed. Qualitative selection of sub-regions was based 
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Fig. 2. MODIS Rrs(645) (derived using the methodology developed in Section 2.1.5) in our four river plume regions of interest during heavy outflow (a) Amazon River, b) Chesa- 
peake Bay, c) Mississippi River, 4) Yangtze River). Sub-regions include the plume waters (black box), estuarine/near-shore/coastal waters (red box), and shelf waters (green 
box). Virtual transects designed for GSD analysis are shown in gray with stations (points) denoted. Each plot shows the date and time of MODIS acquisition and scale. White 
areas over water are masked for clouds and stray light around clouds and land. 
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on visual inspection of satellite imagery for 645). Specifically, sea- 
ward boundaries of typical plume fronts were used to define the 
plume, while the coastal sub-regions were situated in areas found 
to dynamically shift between turbid and clear waters, and continental 
shelf waters were comparatively quiescent and dark at 645 nm. 

2.1.1. Amazon River 

The Amazon River drains a region of South America over 
6 x 10 6 km 2 and has the highest discharge of any river in the world, 
with flow rates between 7.5 x 10 4 and 3.8 x 10 5 m 3 s _l and an 
average discharge of 2.09 x 10 5 m 3 s” 1 (Guyot, Callede, Molinier, 
Guimaraes, & Oliveira, 1997; Moreira-Turcq, Seyler, Guyot, & Etcheber, 
2003) with which it delivers -3.27 x 10 13 g C yr” 1 into the Atlantic 
Ocean. TSM is high throughout the Amazon River network, with some 
regions reaching values over 300 mg l” 1 (Moreira-Turcq et al., 2003). 
Flow is seasonally tied to rain fall, reaching average peak flow in May, 
and nadir flow in November (Zeng, 1999), and is strongly linked to El 
Nino-Southern Oscillation (ENSO) (Molinier et al., 2009). The plume 
of the Amazon River extends several hundred kilometers seaward, 
and— like many tropical regions— is characterized by frequent, heavy 
cloud cover. Hu, Montgomery, Schmitt, and Muller-Karger (2004) stud- 
ied the plume further from shore between 1997-2002 using concurrent 
SeaWiFS satellite ocean color and in situ observations, and found that 
patches of fresher, CDOM-dominated waters originating from the Ama- 
zon and Orinoco extended as far as 2000 km to the north and west of 
the Amazon mouth. Large (400-500 km) eddies were observed moving 
northward from the plume along the coast at a frequency of about 
1-2 months in all seasons. The North Brazil Current and the North Equa- 
torial Counter Current carried colored waters from the plume >3000 km 
into the tropical Atlantic in 1998. Plume areal extent between 2000 and 
2004 mapped by correlating the SeaWiFS-derived detrital and dissolved 
absorption coefficient to sea-surface salinity ranged from 268 x 10 3 km 2 
to 1506 x 10 3 km 2 (Molleri, Novo, & Kampel, 2010). 

713 MODIS Aqua and Terra scenes covering the Amazon plume 
were selected during May and November 2008 (La Nina; peak 
rains), 2009, 2010 (El Nino; less rain), and 2011 (La Nina; http:// 
www.cpc.noaa.gov/products/analysis_monitoring/lanina/ 
enso_evolution-status-fcsts-web.pdf). 

2.1.2. Chesapeake Bay 

Located in the eastern United States and draining an area of about 
1.65 x 10 5 km 2 , the Chesapeake Bay receives the majority of its fresh 
water and nutrient loading from the Susquehanna River (average dis- 
charge 1.06 x 10 3 m 3 s” 1 (Gellis, Banks, Langland, & Martucci, 2005)) 
at the head of the Bay (Harding, Itsweire, & Esaias, 1994). Sediment 
loading combined with high concentrations of CDOM from the Susque- 
hanna and other tributaries lead to diffuse attenuation of light in the 
water column, I< d (490), as high as -10 m” 1 (Wang, Son, & Harding, 
2009), which in turn degrades benthic habitat critical to local fisheries. 
TSM concentrations of 10-50 mg l” 1 are not uncommon at the head of 
the Bay, but are generally <10 mg l” 1 in the central and lower Bay 
(Ondmsek et al., 2012). 

The use of a specialized atmospheric correction for turbid waters 
using both SWIR and N1R bands for aerosol correction has yielded en- 
couraging results in the Chesapeake Bay for retrieving l< d (490) (Wang 
et al., 2009). However, an evaluation of the method for SWIR aerosol 
correction as implemented in SeaDAS software was criticized in 
(Werdell, Franz, & Bailey, 2010) for yielding a much larger range of 
variability and higher frequency of negative water leaving radiance 
retrievals; results consistent with increased SWIR atmospheric 
correction noise as discussed in Section 1.3. 

630 MODIS images were selected from average peak (March) and 
nadir (September) flow periods between 2008 and 2011 based on Unit- 
ed States Geological Survey (USGS) river gauge data collected at the 
Conowingo Dam (http://waterdata.usgs.gov/usa/nwis/uv701578310). 


2.1.3. Mississippi River 

The Mississippi River drains 41% of the continental United States 
(3.27 x 10 6 km 2 ) with an annual discharge rate of 1.84 x 10 4 m 3 s” 1 
and averages a flow rate of 1.35 ± 0.2 x 10 4 m 3 s” 1 (Milliman & 
Meade, 1983). Plume waters measured from satellite have been 
shown to reach as far as the Gulf Stream (Hu et al., 2005), and the 
plume extent (TSM > 5 mg l” 1 ) following a flooding event in 2008 
was estimated from satellite to be as high as 5859 km 2 , or about double 
that of the climatological average (2002-2008) (Shi & Wang, 2009b). 
The estimated 2.1 x 10 s tons (Milliman & Meade, 1983) of sediment 
delivered by this flow to the shelf annually are accompanied by high in- 
organic nutrient loading contributing to recurrent hypoxia in the bot- 
tom waters of the shelf (Rabalais et al., 1996). Average peak and nadir 
flow periods for our study period were May and September, respective- 
ly, based on USGS gauge data (http://waterdata.usgs.gov/usa/nwis/uv? 
site_no=07374000). 803 MODIS images were selected from the Missis- 
sippi region in the period between 2008 and 2011. 

2.1.4. Yangtze River 

The Yangtze River drains about 1.94 x 10 6 km 2 and has an annual 
discharge of 2.85 x 10 4 m 3 s” 1 (Milliman & Meade, 1983) making it 
the fourth largest in the world. While suspended sediment delivery to 
the shelf (-5.0 x 10 8 tons annually) is principally controlled by river 
flow (average 3.0 x 10 5 m 3 s” 1 ), spring and neap tides play a significant 
role (Milliman, Huang-ting, Zuo-sheng, & Mead, 1985). Seasonal river 
flows and sediment discharges peak in July with the Southeast Asian 
summer monsoon, and are generally near nadir in January (Beardsley, 
Limeburner, Yu, & Cannon, 1985) when prevailing winds turn offshore. 
Surface turbidity as measured by I< d (490) estimates from MODIS have a 
nearly opposite pattern on the nearby shallow shelf owing to high winter 
winds and a breakdown in stratification leading to resuspended particles 
mixing to the surface (Shi & Wang, 2010). 

The plume of the Yangtze is defined here to include the Hangzhou 
Bay at the mouth of river, as well as adjacent shallow (-20 m) shore- 
line regions where turbidity is driven by a combination of these fac- 
tors (Fig. 2). 757 MODIS images were selected from the region 
during January and July 2008-2011. 

2.1.5. Atmospheric correction and quality assurance 

SeaDAS software (http://seadas.gsfc.nasa.gov/) provides a reliable, 
user-friendly interface for atmospherically correcting ocean color data 
from several instruments over the vast majority of the world's oceans. 
However, several challenges must be overcome in order to accurately 
process MODIS ocean color imagery in highly turbid zones— the most sig- 
nificant being the atmospheric correction. Traditionally, visible (VIS, 
400 nm-700 nm) bands are atmospherically corrected using two bands 
in the NIR (748 nm (band 15, 10 nm bandwidth) and 869 nm (band 
16, 15 nm bandwidth) for MODIS) to remove contributions to the TOA 
reflectances attributable to aerosols. Standard atmospheric correction al- 
gorithms for the global oceans assume negligible water leaving radiance 
in the NIR due to the strong absorption coefficient for water in this region 
of the spectrum (the “black pixel” assumption (Gordon & Wang, 1994)). 
More recently, this approach has been refined in SeaDAS to include an it- 
erative correction with a bio-optical model for estimating L^fNIR ) (Bailey 
et al., 2010). However, in highly turbid coastal waters, strong particulate 
backscattering often leads to relatively high l^fNIR), thereby violating 
the NIR black pixel assumption, and often leading to non-convergence 
of the iterative bio-optical optimization in extreme conditions, as we 
show below. Furthermore, NIR bands on MODIS sensors Aqua and Terra 
both saturate at TOA radiances above -2.8-3.45 mW cm” 2 pm” 1 sr” 1 
(748 nm) and -1.9-2.45 mW cm -2 pm” 1 sr” 1 (869 nm), which 
leads to data dropouts over the most reflective plume regions 
(>-30-35 mg l” 1 TSM). 

The default cloud mask in SeaDAS also depends on cloud albedo in 
the 869 nm band, and even in areas where it is not saturated, the 
default threshold for cloud identification (cloud albedo >0.027) is 
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subject to significant error due to non-negligible water leaving radi- 
ance at 869 nm in turbid waters (Wang & Shi, 2006). 

A method for avoiding violation of the NIR black pixel assumption 
has been developed which utilizes MODIS bands in the SW1R 
(1240 nm (band 5, 20 nm bandwidth) and 2130 nm (band 7, 50 nm 
bandwidth) in MODIS) (Wang, Son, & Shi, 2010). At the longer of 
these two wavelengths, water leaving radiance is effectively zero even 
in highly turbid waters because the absorption coefficient for water in 
the SW1R is more than two orders of magnitude higher than in the 
NIR (Shi & Wang, 2009a). Unfortunately, the standard application of 
this atmospheric correction in SeaDAS requires the use of the NIR 
bands (as well as additional requirements on Chi and nL w ) to establish 
a turbidity index above which the SWIR correction is applied, and in 
highly turbid waters such as those often found in our areas of interest, 
these bands are saturated. To further complicate matters, the 
1240 nm channels in both the Aqua and Terra instruments— originally 
designed for land applications— have low SNR as well as malfunctioning 
detectors, leading to increased noise and missing data in the final, atmo- 
spherically corrected L2 product. The 2130 nm channel on Terra is also 
subject to data dropouts due to an instrument crosstalk issue. 

Other challenges to processing MODIS imageiy in highly turbid wa- 
ters include the default masks in SeaDAS L2 file generation scripting 
(the 12gen module) for high light and stray light. Stray light is light 
reflected into the detector from nearby bright land or cloud pixels. 
These masks often hide proportionately large and sometimes critical 
areas of nearshore, estuarine, lake, and inland waters. High light 
masking is based on saturation or near-saturation in the MODIS NIR 
bands, and therefore cannot be used in highly turbid environments 
(i.e. they may saturate before other sensors/bands of interest). No 
stray light correction is currently attempted in SeaDAS for MODIS. 
Stray light masking— while not active by default— is nevertheless impor- 
tant for eliminating light contamination in pixels adjacent to bright tar- 
gets. When activated, the default mask of 5 along-track pixels by 7 
across-track pixels around bright land and cloud/ice pixels can signifi- 
cantly reduce data coverage, and may be impractical in nearshore, estu- 
arine, and inland waters. 


Experimentation with various atmospheric correction and masking 
permutations guided us in our approach to MODIS L2 file generation 
in SeaDAS for our river plume regions of interest. R rs products from 
each method were evaluated for signal contamination (e.g. excessive 
stray light and other artifacts resulting from inadequate masking), 
poor atmospheric correction (e.g. high frequency of negative or other- 
wise unrealistic reflectances), and poor retrieval (e.g. discontinuities 
and values outside of the range of field measurements found in the 
literature). It should be noted that no comprehensive, quantitative 
validation of R rs retrievals is attempted here. We insure a continuous 
transition in each scene between products derived using standard 
methods in darker offshore waters (as described above and in (Bailey 
et al., 2010; Gordon & Wang, 1994)) and those derived using methods 
optimized for turbid pixels (described below). Furthermore, we com- 
pare retrievals of TSM using our methodology to field measurements 
for one particular example in the Chesapeake Bay following major 
flooding caused by Hurricane Irene and Tropical Storm Lee in 2011. 
The primary goal of this study— evaluating optimal GSD in ocean color 
imagery— depends principally on the interpixel variability in ocean 
color products rather than their absolute magnitude. Nevertheless, 
qualitative and quantitative steps are described below which insure 
that retrievals of R rs are within reasonable ranges. 

The first processing step ( see Fig. 3 for a processing flow-chart) was to 
estimate R^(645) for each scene with a methodology optimized for high 
turbidity pixels (hereafter referred to as HT). Unix shell batch scripts 
called the SeaDAS 6.4 subroutine 12gen, forcing the use the SWIR aerosol 
correction (i.e. parameter options aer_opt = —1, aer_wave_short = 
1240, aer_wave_long = 2130). To mitigate data loss due to bad detectors 
and crosstalk issues in the SWIR bands, and to reduce the impact of low 
SNR in these bands, the 1240 and 2130 bands were smoothed by a 
3x3 pixel filter prior to atmospheric correction. This is accomplished 
by defining a custom msll2_filter.dat parameter file (found in the SeaDAS 
SOCDATAROOT directory) and adding lines for bands 14 and 16 (i.e. 
ltmean, 14, 3, 3, 1; ltmean 16, 3, 3, 1). High light masks were disabled, 
and stray light masks were set to a 3 x 3 array around land and clouds 
(i.e. maskland = 1, maskhilt = 0, maskstlight = 1, maskcloud = 1, 
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Fig. 3. Processing flowchart for deriving a blended satellite product using a method optimized for high turbidity (HT) pixels and SWIR aerosol model selection as well as low tur- 
bidity pixels (LT) and bio-optical iteration of the NIR aerosol model selection as described in Section 2.1.5. Shaded regions indicate software used during process flow (i.e. SeaDAS 
and Matlab). 
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filter_file = $0CDATAR00T/modisa/msll2_filter_custom.dat in the 
scripted call to 12gen, and stlight, 10, 3, 3, 1 in the custom filter file). The 
2130 nm SWIR band was used to mask clouds with a threshold albedo 
of >0.018 (cloud_thresh = 0.018). Unrealistic R,s(645) retrievals 
(Rrs( 645) < = —0.02 sr _1 or R rs (645) > 0.12 sr~ ') were discarded. 
Scenes with >85% cloud cover were also discarded as they have the po- 
tential to introduce excessive noise due to cloud shadow and error in 
cloud masking. The remaining scenes were output at QKM, HKM, and 
1 KM resolution, the latter two of which are derived from native resolution 
in SeaDAS using pixel aggregation and averaging (http://mcst.gsfc.nasa. 
gov/sites/mcst.gsfc/files/file_attachments/Ml 054.pdf). 

For darker waters, a variation on the standard OBPG NIR atmo- 
spheric correction with bio-optical iteration was considered more ac- 
curate due to the low SNR in MODIS SWIR bands (see Section 1.2) 
(Werdell et al., 2010). As in HT, this low-turbidity approach (LT) 
uses 2130 nm for cloud detection, masks stray light but not high 
light, and removes outliers based on the reflectance thresholds 
given above. Histogram analysis of retrievals for R rs (645) generally 
showed a population of slightly negative pixels in the darkest regions 
of the image, indicating a slight systematic bias due to sensor 
radiometric calibration or overcorrection for aerosols. The median of 
negative retrievals (not less than — 0.02 sr _1 ) was therefore calculat- 
ed for each scene (generally <sc— 0.001 sr -1 , but occasionally as high 
as —0.002 sr" 1 ) and subtracted from all pixels processed in LT. This 
slight offset— used to compensate for minor calibration and atmo- 
spheric correction error in SeaDAS— has a negligible effect on turbid 
waters, but it improves coverage offshore because negative retrievals 
are subsequently eliminated. 

Next the two 12gen schemes (LT & HT) were merged. Histogram 
comparisons of low reflectance pixels in HT with those in LT revealed 
small biases in the SWIR approach (HT), possibly owing to the lack of vi- 
carious calibration in MODIS SWIR bands. An offset was calculated to 
match the modes of HT and LT for low reflectance pixels (i.e. those 
not flagged as HILT) and added to the HT scenes. These offsets were 
found to be very small generally ( — 0.002 to 0.002 sr _1 ) and more fre- 
quently slightly negative, but helped insure a continuous geographic 
transition between algorithms. Finally, scenes were merged based on 
a threshold R rs (645) value of 0.01 sr -1 (corresponding approximately 
to TSM of 9.5 mg l -1 (Miller & McKee, 2004)) resulting in a “L2B” file 
(Fig. 3); i.e. LT values were applied to all pixels below this threshold, 
and HT values were retained elsewhere. It is noteworthy that threshold 
selection is applied to HT imagery— replacing low turbidity pixels with 
LT values— rather the other way around, which would be subject to 
large errors for regions with bio-optical iteration failure. Bio-optical 
iteration failure occurs when the estimation of NIR water leaving 
radiance for selecting the appropriate aerosol model fails to converge 
as described in (Bailey et al., 2010), and results in an ‘atmospheric cor- 
rection warning’ flag and— in our experience in highly turbid waters— 
significant underestimates of remote sensing reflectance in the red, as 
shown in Fig. 4. 

Prior to merging HT and LT scenes, imagery was rejected in which 
>85% of ocean pixels were masked for cloud. While this eliminated 
most heavily clouded scenes potentially susceptible to higher interpixel 
error in ocean color due to cloud and stray light effects, in many 
instances it failed to eliminate scenes with poor coverage of our 
sub-regions of interest and transects resulting from only partial granule 
coverage of the region ( i.e. cloud cover was evaluated as a percentage of 
available non-land masked pixels in an image subset to our region, not 
of potential pixels covering the entire region). All scenes were therefore 
visually inspected and discarded for poor coverage deriving from tran- 
sects and regions of interest lying off-granule (the overwhelming pro- 
portion of manually rejected scenes), and for heavy, patchy cumulous 
cloud cover, and— in a few instances— unrealistic mode offsets deriving 
from too few LT pixels available in the scene. Specifically, 115 of the 
remaining 700 scenes were discarded from the Aqua data set, and 
148/757 from Terra. 


2.2. Finding optimal GSD 

The procedure for finding optimal GSD followed the theoretical ap- 
proach outlined in Section 1. Tabulated data for SNR as a function of 
TOA radiances for each MODIS Aqua band and detector were provided 
for this study by the NASA MODIS Characterization Support Team 
(MCST) and were applied to Eq. (5). Their methodology for deriving 
on-orbit SNR is described in (Xiong, Sun, Barnes, & Salomonson, 
2010). At the time of writing, tabulation was not available for Terra, 
but based on preliminary results, it does not appear that differences 
will be high for band 1 , and the Aqua values are applied here to Terra 
as well. Results showed small differences between the two instruments, 
with higher interpixel variability in Terra (e.g. Section 3.2), which may 
be a result of an overestimation of SNR. 

For each QKM scene in each region, the transects shown in Fig. 2 were 
traversed using software developed in Matlab (www.mathworks.com). 
At each station, aggregations of 2 x 2,3 x 3,4 x 4, etc., pixels were eval- 
uated until the coefficient of variability in Rrs(645), oj, was greater than 
the noise, a t . For each image, the actual resolution was calculated, and 
any scenes with pixels larger than 500 m were rejected, as was any sta- 
tion with CTj > 1.5 or where >50% of pixels were undefined due to 
masking for clouds, land, etc. a; > 1 .5 is defined as the standard deviation 
of the array equal to 1 50% of its average value, which is only likely if some 
of the pixels are unflagged land, cloud, or straylight Once the interpixel 
variability threshold was met, the optimal GSD was defined as the aver- 
age between the size of the inconclusive array (Oj < a t ) and the size of 
the array which marked the upper limit of a GSD which could resolve sig- 
nificant differences in R re ( 645) (Oj > a t ). 

2.3. Semivariance 

Semivariance analysis was conducted only within the plume regions 
shown in Fig. 2 and in the Sargasso Sea. Analysis relied upon Matlab code 
adapted from W. Schwanghart and available at Mathworks® (http:// 
www.mathworks.com/matlabcentral/fileexchange/20355), which test- 
ed 5000 random pairs of pixels in any orientation over a range of lags 
from 1 pixel distance to the maximum distance available in the subre- 
gion. Experimental variograms were then fitted using least square differ- 
ence minimization to the spherical model: 

y(h) = c 0 + Cl (l .5 -0.5 Qj 3 V (6) 

where c 0 , Ci, h, and a are the nugget, sill, distance, and range, respectively 
(see Section 1.2). Every experimental semivariogram was visually 
inspected, and the range and sill of the smallest patch was manually se- 
lected based on the flattening of the experimental semivariance at 
shortest h prior to fitting. Quality of fits of experimental data to the 
spherical model were evaluated qualitatively and rejected in cases 
where the modeled nugget was significantly different from the experi- 
mental, and where the shape could not be reproduced with a spherical 
model. Spherical fits to the experimental data are used here only to aid 
in visualization of results, and key parameters such as nugget, sill, and 
range are all derived from experimental results. 

3. Results and discussion 

3.1. Tropical Storm Lee 

Two major tropical systems impacted the Chesapeake Bay region 
in late summer 2011 providing a good opportunity within the scope 
of this study to demonstrate some of the benefits of the methodology 
for turbid waters developed here. The first was Irene, which passed 
over the southern Mid-Atlantic States on August 27-28 as a Category 
1 hurricane, bringing as much as 30 cm of rain to some portions of 
Maryland (http://www.nhc.noaa.gov/data/tcr/AL092011_Irene.pdf) 
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Fig. 4. An example image from the upper Chesapeake Bay on September 11, 201 1 following Hurricane Irene and Tropical Storm Lee (MOD1S granule A201 12541 80500) comparing 
various standard N1R and SW1R-N1R approaches to atmospheric correction in SeaDAS to the method used here, a) an enhanced true color image (HKM), b) standard bio-optical, 
iterative N1R process in SeaDAS (i.e. aer_opt = — 3) with default masks (land, cloud/ice, and high light (HILT) at QKM, c) same as b with HILT and erroneous cloud masks removed 
to show underlying error in the bio-optical iteration within the plume, as well as atmospheric correction failure (ATMFA1L) owing most likely to N1R band saturation in extreme 
turbidity, d) standard SWIR-NIR product (i.e. aer_opt = — 9) showing nearly identical results to c— though noisier in dark waters— due to failure to switch to SWIR bands for aero- 
sol correction within the plume. The failure to switch probably arises from erroneously low nL w (see c). e) Results from this study processed as described in Section 2.1.5 using SWIR 
bands in the plume, N1R bio-optical iteration in darker waters, and custom cloud, HILT, and stray light masking throughout with no reliance on N1R bands. 


and saturating the surrounding watershed. While increased turbidity 
was detectable in MODIS imagery in the days following Irene due to 
the increased run-off and resuspension, the most dramatic impacts 
came over a week later when the remnants of Tropical Storm (T.S.) 
Lee deposited an addition 53 cm of rain over parts of the region 
based on National Oceanic and Atmospheric Administration (NOAA) 
gauge data (http://www.nhc.noaa.gov/data/tcr/AL132011_Lee.pdf). 
Flooding reached 6.7 m above flood stage on the Susquehanna 


River, and on September 6, the Conowingo Dam at the head of Ches- 
apeake Bay was opened to relieve the flooding, increasing flow from 
-2.3 m 3 s -1 to 17.0 m 3 s -1 (http://waterdata.usgs.gov/usa/nwis/ 
uv?01578310) and releasing over the next several days a turbidity 
plume extending nearly 150 km down the Bay (Figs. 4 and 5). 

On September 12 and 13, field measurements of TSM were made 
within the plume (Fig. 5) following the methods described in 
(Ondrusek et al., 2012). Applying the algorithm of (Ondrusek et al., 
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Fig. 5. Upper panels show field collection sites for September 12-13, 2011 overlaid on Aqua retrievals of R rs (6 45) derived as per Section 2.1.5. White pixels represent masking for 
clouds and stray light around clouds and land. Image days, times, and resolution are noted in the northwestern corners of each. Match-ups between TSM retrieved using (Ondrusek 
et al., 2012) and field measurements are shown in the bottom plot together with 1 : 1 line (dashed), fit line (solid), and regression statistics including r 2 , bias, root mean square error 
(RMSE), percent difference (PD), slope and intercept of the fit. Retrievals from Aqua are in red, Terra are black. 


2012), TSM was retrieved from MODIS imagery processed as 
described in Section 2.1.5 and matched to field stations based on a 
3x3 pixel array within 3.5 h. Given the cubic nature of the TSM algo- 
rithm, and the fact that it was tuned to considerably lower values of 
TSM than those shown here, the retrievals in the most reflective waters 
tended to be overestimated. For this reason retrievals of >150 mg l -1 
(two stations) were excluded. The lower panel in Fig. 5 shows reason- 
able agreement between retrievals and field measurements, although 
the sample size is quite limited, and correlations may be hampered 
by the lack of any data in the low turbidity portion of the graph 
(<~50 mg L -1 ). 

Figs. 6a and b show the results of GSD and semivariogram analysis, 
respectively. Within the plume, the size of the pixel arrays around vir- 
tual stations showing interpixel variability above the noise threshold 
(i.e. variability in ocean color) is never more than 5x5, and most 
often just 2x2 (nominally 250 m pixels). An increase in necessary 
array size (decrease in resolution) can be seen in the estuary (here 
the central and lower Chesapeake Bay), particularly in image 
A2011255. On the shelf, resolutions as low as 4 km were sometimes 
adequate (16x16 pixels), although this varied highly from image to 


image. Within the plume itself, the smallest patches had ranges of 
-12-20 km (Fig. 6b), and the semivariogram nuggets were <~0.15, 
indicating that >85% of the interpixel variability within the patch 
was detectable at QKM. 

3.2. Optimal ground space distance 

Optimal GSD results from virtual transects in each region are 
presented in Fig. 7 and Table 1 . In general, the pixel size required to 
capture variability in ocean color increases from inshore to offshore 
as expected. GSD requirements have large spatial and seasonal vari- 
ability among the over 19,000 virtual stations remaining after quality 
in the imagery was assured— and often are limited by the 375 m 
lower bound of our analysis— but average between 522 ± 57 m in 
combined plume regions and 1354 ± 92 m in the open ocean. 

For nearly all regions, the results for the sub-regions we selected (i.e. 
plume, nearshore, and shelf waters) were statistically different 
(ANOVA, p <sc 0.01, see Table 1). An exception was found between 
nearshore and shelf waters in the Amazon region. Within sub-regions, 
some significant differences (eight of 12) were found between Aqua 
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a Chesapeake GSD Transect; T.S. Lee 



b Chesapeake (QKM); T.S. Lee 



Fig. 6. Examples of both optimal GSD results (a) and semivariograms (b) from the 
Chesapeake Bay region following the passage of Hurricane Irene and Tropical Storm 
Lee. Colors in both a and b represent individual MODIS QKM images (legend in b). Col- 
ored lines in a are linear fits to the data from all virtual stations (colored dots) along 
the transect showing general trends toward increasing GSD with distance from the 
plume. Each dot in a represents the pixel array size at which the interpixel variability 
in ocean color exceeds that of instrumental or atmospheric noise. Coarser resolution 
would overlook these differences in ocean color. In b, experimental semivariograms 
(dashed) are overlaid on spherical fits (solid). Nuggets all fall below -0.15 indicating 
that QKM imagery is capturing >85% of the overall variability within the patch size 
(i.e. range, -12-20 km). 

and Terra satellite results, with Terra falling 15% lower on average than 
Aqua. The discrepancy likely derives from higher noise levels in the 
Terra instrument which are not accounted for in the on-orbit noise es- 
timation data available (see Section 2.2). 

Only four of the 12 comparisons between GSD during peak and 
nadir river flow periods showed significant differences: in the 
plume sub-regions of the Chesapeake and Mississippi, and the near- 
shore and shelf sub-regions of the Yangtze. 

A dip in the optimal GSD and reduction in inter-image variability 
is visible in the plume region of the Amazon between about 200 
and 350 km (Fig. 7) indicating a stronger and less transient gradient 
in ocean color. This suggests that the geographic and temporal extent 
of this plume is comparatively stable relative to the other regions 
studied. 


It is not immediately apparent what drives the spatial variability we 
see in the Sargasso Sea, where optimal GSD varies between station av- 
erages of -750-4400 m. The region was considerably larger than the 
others (10° latitude by 10° longitude) and prone to heavy cloud cover. 
Of the over 1100 granules with at least some coverage of the region dur- 
ing the periods chosen, only 109 were retained for analysis after elimi- 
nation for 1 ) having no open sky pixels over the virtual transect— which 
stretched the extent of the region from the southwest to the northeast— 
or 2) for excessive cloud cover (>65% of pixels). Remaining images were 
never completely cloud free (12% minimum), and averaged 46% 
clouded overall. If stray light masking around clouds was insufficient 
and/or cloud shadows a factor, this may have been mistaken by the 
GSD analysis for interpixel variability in ocean color, thereby reducing 
the GSD results in cloudy scenes. 

3.3. Semivariogram analysis 

Semivariograms— measures of the variability and spatial extent of 
ocean color patches— differed significantly in the four river plume re- 
gions generally as a function of the relative size of the rivers themselves 
(Table 2, Figs. 8-10). The largest plumes at the Amazon and Yangtze 
Rivers had the largest coherent patches (87 ± 21 km, 62 ± 4 km, re- 
spectively), followed by the Mississippi (43 ± 4 km) and the Chesa- 
peake (19 ± 6 km). Due to the large patch sizes, semivariograms 
were largely insensitive to the resolution used to measure them. 
These patches do not represent homogeneous waters, but spatial scales 
beyond which variability in ocean color does not increase— at least for 
some distance. Variability within patches also depended on river size 
with sill values declining from 1.25 ± 0.34 sr -1 in the Amazon 
plume, to 0.61 ± 0.30 sr -1 at the Yangtze, and 0.32 ± 0.26 sr -1 and 
0.12 ± 0.09 sr" 1 at the Mississippi and Chesapeake, respectively. By 
contrast, the proportion of variability resolvable at a given resolution 
(i.e. 1— nugget), varied from 91% at the Chesapeake to 95%, 96%, and 
98% at the Yangtze, Mississippi, and Amazon, respectively. The nugget 
was the only parameter that showed significant differences between 
seasons, with more variability going undetected during nadir flow sea- 
sons, particularly in the Chesapeake (Table 2). In the Chesapeake, large 
nuggets during nadir flow of the Susquehanna River may result from 
the influence of the many other un-gauged tributaries feeding into the 
Susquehanna plume. The spherical model (Eq. (6); used only for data 
visualization) fit the empirical data about 70% of the time based on 
the criteria outlined in Section 2.3, and was more likely to fit well if 
the range was short. 

Results from our control region in the Sargasso Sea were consider- 
ably different from those in the plume regions. Patch size was much 
larger and more variable at 180 ± 77 km, while the variability within 
the patches was 2-3 orders of magnitude lower (4.1 x 10“ 7 ± 
5.0 x 10“ 7 sr" 1 ) than in the plumes, and degree of unresolved vari- 
ability increased to 34% on average. Unlike other regions, the sill was 
highly sensitive to resolution (Table 2), although the reasons for this 
are not immediately apparent, and it was not reflected in other 
semivariance parameters. 

Figs. 8-10 show experimental and fitted semivariograms in greater 
detail for each region at HKM (chosen for being closest to GSD results 
(Table 1), and not significantly different from other resolutions 
(Table 2)) together with distributions of the parameters summarized 
in Table 2. Due to the similarities in semivariograms discussed above, 
these figures are separated into the largest rivers (Fig. 8), moderately 
large rivers (Fig. 9), and offshore (Fig. 10). Overall, the smaller patches 
(lower range) tended to fit better to the spherical model (Eq. (6); 
dashed green and brown lines). Experimental semivariograms (gray 
lines) were considerably noisier in the Sargasso Sea owing to the vari- 
ability in the reflectance signal being extremely low, but still matched 
the spherical model in the majority of cases. The low sill values found 
offshore were also likely the reason for the notably higher degree of 
unresolved variability as a proportion of the sill. Had the nuggets not 
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Amazon Transect Chesapeake Transect 



Sargasso Transect 



Fig. 7. Optimal ground space distance for each region as a function of distance along the river plume transects in Fig. 3 and the Sargasso Sea (note changes in scales of GSD and distance). 
Sub-regions are summarized in the box-whisker diagrams to the right of each transect. As with all box-whisker diagrams in this manuscript, the thin central line shows the population 
median, boxes extend to the first and third quartiles, vertical lines extend from boxes to -99% of the data, and extreme outliers are shown as red plus symbols (partially truncated to pre- 
serve scale on occasion). The thick, solid red lines along the transects represent the median values of all satellite images for the region over all seasons ( see text for numbers of images), and 
the shaded regions are one standard deviation. 


been normalized to sill values in the offshore region, their magnitudes 
would have also been orders of magnitude lower than those found in 
the river plumes. 

It is perhaps somewhat surprising that neither satellite resolution 
nor season generally play a significant role in the spatial scales of 
patchiness in the plume regions. In part, this is attributable to the 
higher degree of variability among all images as shown in the 
box- whisker plots in Figs. 8-10. Another factor— at least as far as res- 
olution comparisons are concerned— may be the size of the patches 
themselves. Patches were identified visually as the smallest ranges 
within which the semivariance plateaued with a discernible sill, and 
yet they were always much larger than many dozens of pixels; the 
smallest average for region, season, and resolution being 12 km at 
QKM, or -48 pixels. This likely explains why results did not vary no- 
tably between 1KM, HKM, and QKM. 

4. Summary 

Improving spatial resolution in future ocean color satellite sen- 
sors such as GEO-CAPE incurs additional expense and complexity of 


instrument design, and should be justified by increased capabilities. 
Most of the data available from low Earth orbiters has historically 
been collected at - 1 km resolution, although exceptions are available 
from MERIS (300 m), some bands of MODIS (250 m and 500 m), and 
others. Previous studies have looked at high resolution space- and 
air-borne imagery (generally a small number of scenes) to understand 
what resolution is adequate for mapping biogeochemical property dis- 
tributions in coastal and inland waters. Here, we have expanded upon 
and developed analytical methodologies for satellite ocean color pro- 
cessing and for determining optimal GSD, applying them to four years 
of MODIS imagery (over 3000 images) from the red 645 nm band at 
QKM, HKM, and 1KM resolution. 

For the purpose of demonstrating the method used for satellite 
imagery processing as well as the spatial statistics tools employed 
here, the period following two tropical systems crossing the Chesa- 
peake Bay was examined in detail. Next, four river plumes at the Am- 
azon, Chesapeake, Mississippi, and Yangtze Rivers were studied with 
adjacent waters, and an offshore region was shown for comparison. 
Regions were sectioned into sub-regional environments, transected 
by virtual stations used to extract satellite imagery, and analyzed for 
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Table 1 

Median optimal GSD for each region as a function of season, satellite, and sub-region. 
One-way ANOVA tests between adjacent values indicate the likelihood that popula- 
tions have a significantly different (red cells) or similar (green cells) median. Compar- 
isons with p-values > 0.05, between 0.01-0.05, 0.001-0.01, and <0.001 are labeled NS, 
>, <, and <<, respectively. All comparisons between Plume and Shelf waters were signif- 
icantly different (p <c 0.01). The summary at the bottom of the table includes the 
mean GSD and numbers of virtual stations (N) used in each subregion. 


Plume <-> Nearshore <-> Shelf Offshore 
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Aqua 

i 
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1262 

NS 

1446 

Summary 


522±57 

630±43 

753±106 

1354 ±92 

N 


6226 

9299 

3232 

471 


optimal ground space distance, spatial scales of variability, and the 
proportion of variability detectable within plume patches. 

Unsurprisingly, the plumes themselves were found to be the most 
spatially heterogeneous with scales of patchiness (i.e. ranges) on the 
order of -20-90 km compared with -180 km offshore. In plumes, 
these scales correspond to processes such as fronts, eddies and fila- 
ments, synoptic storm outflows and resuspension, and internal tides, 
while in coastal waters they may also correspond with phytoplankton 
blooms, inertial, internal and solitary waves. Patch sizes found offshore 
correspond to the scale of larger phenomena such as surface tides and 
seasonal changes in mixed layer depth and biomass cycles (Dickey et 
al., 2006). Patches were large enough that the resolution (QKM, HKM, 
or 1KM) used to calculate the semivariance function was not a signifi- 
cant factor. The ranges in the river plumes were roughly proportional 
to the sizes of the rivers studied, as were the degrees of variability (i.e. 
sills) within the patches. The proportions of the variability resolved 
(i.e. nuggets) in the patches were generally high at all resolutions mea- 
sured, and varied between plume regions from 85% to 99%, with the 
lower values found in the Chesapeake Bay. 

Significant differences exist between the spatial resolution required 
to detect interpixel ocean color variability between the three coastal en- 
vironments (-520 m, 630 m, 750 m, and at the plumes, nearshore, and 
shelf, respectively) when both the peak plume and nadir plumes sea- 
sons are taken together. Separately, the nadir river flow seasons tend 
to require slightly lower GSD ( 7%-l 2% in the plume, mean across subre- 
gions 4%) than peak flow seasons. Results from Terra were generally 
somewhat lower than those from Aqua (-12% lower GSD in the 
plumes) possibly as a result of higher, unaccounted for noise in the 
instrument. 


Table 2 

Experimental semivariance component data for each region, season, and resolution. Few 
significant differences were found in the between Aqua and Terra or between 1KM, 
HKM, and QKM within a region (p > 0.01). Peak nugget populations tended to be slightly 
lower than nadir (ll,p« 0.01), except in the Chesapeake where the difference was pro- 
nounced (10%). N is the number of semivariograms used in each region/season/resolution 
and n is the number successfully fit to the spherical model. 



Season 

Nugget 

[%] 

Range 

[km] 

Sill 

[xl0' 4 l/sr] 

N(n) 

Amazon 

1KM 

Peak 

2 

110 

1.630 

34(16) 



Nadir 

1 

50 

0.692 

51(33) 


HKM 

Peak 

2 

93 

1.275 

36(22) 



Nadir 

1 

79 

1.086 

69(47) 


QKM 

Peak 

2 

90 

1.271 

54(41) 



Nadir 

3 

98 

1.537 

90(58) 

Chesapeake 

1KM 

Peak 

3 

26 

0.187 

55(26) 



Nadir 

15 

16 

0.039 

47(41) 


HKM 

Peak 

4 

23 

0.216 

73(45) 



Nadir 

14 

13 

0.041 

35(34) 


QKM 

Peak 

5 

23 

0.216 

75(29) 



Nadir 

13 

12 

0.045 

37(14) 

Mississippi 

1KM 

Peak 

2 

44 

0.518 

136(91) 



Nadir 

4 

51 

0.067 

123(82) 


HKM 

Peak 

3 

41 

0.572 

140(115) 



Nadir 

4 

43 

0.095 

110(87) 


QKM 

Peak 

3 

41 

0.580 

143(124) 



Nadir 

6 

40 

0.103 

126(109) 

Yangtze 

1KM 

Peak 

4 

68 

0.334 

36(28) 



Nadir 

6 

59 

0.888 

82(60) 


HKM 

Peak 

4 

65 

0.329 

39(32) 



Nadir 

4 

59 

0.906 

81(63) 


QKM 

Peak 

4 

65 

0.329 

39(19) 



Nadir 

5 

58 

0.848 

83(47) 

Sargasso 

1KM 

May 

25 

284 

0.00008 

45(27) 



November 

38 

274 

0.00011 

24(17) 


HKM 

May 

44 

151 

0.00160 

45(43) 



November 

32 

136 

0.00270 

20(18) 


QKM 

May 

36 

115 

0.00700 

52(52) 



November 

30 

122 

0.01295 

32(32) 


The offshore region exhibited variability several orders of magnitude 
lower than in the plumes, and that variability was not as well resolved 
(56%-75%). However, a considerably larger pixel size (-1400 m) was 
found to be adequate for resolving interpixel variability. 

The approach suggested here for processing MODIS imagery with 
SeaDAS in highly turbid regions overcomes the problems diagramed in 
Fig. 4 such as over-masking, erroneous cloud detection, iteration failure 
of the aerosol algorithm, and bad switching between SWIR and NIR 
aerosol models, and suggests a future avenue for quantifying pixel- 
specific A Rrs at ocean color bands other than 645 nm in MODIS 
(i.e. Eqs. (1 )-(5)). This approach mitigates problems found using the op- 
erational SeaDAS SWIR-NIR approach in the Chesapeake Bay (Werdell et 
al., 2010) (i.e. higher frequency of negative retrievals and increased noise 
in retrievals) by using a dark water mode offset on SWIR retrievals to ac- 
commodate the lack of vicarious calibration, by unmasking high turbidity 
waters previously hidden by erroneous cloud flags and NIR band satura- 
tion, and by only applying SWIR-selected aerosols in these turbid waters 
where SNR in band 1 MODIS is less of an issue. During the flooding 
following T.S. Lee in the Chesapeake, several clear days provided good 
satellite observing, and field data were collected by associates in the 
resulting plume. Satellite retrievals of TSM showed reasonable agree- 
ment with those collected in the field, though the dynamic range did 
not include low turbidity measurements, the sample size was small, 
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Fig. 8. Semivariograms and related parameters from HKM MODIS imagery at the two larger river plume regions, the Amazon and the Yangtze. The fitted semivariogram curves 
(dashed lines in central area of each plot) are normalized to sill values; non-normalized sill distributions of the experimental data are shown in right-hand box-whisker diagram 
(see Fig. 6 caption for details regarding box-whisker diagrams). Fitted peak (brown dashed lines) and nadir (green dashed lines) river flow periods are overlaid on experimental 
data (grey lines) where they matched well, which happened more often when the range was short (<-100 km). Nugget distributions of non-fitted, experimental data (also nor- 
malized) are shown in the left-hand box-whisker plots by season, and ranges in the bottom box-whisker plots. Differences between seasons and MODIS instruments were small 
(see text). 


and correlations were not strong (r 2 = 0.33). Remote sensing reflec- 
tance retrievals show realistic geographic patterns throughout the 
plumes and darker waters with no obvious discontinuities, little evi- 
dence of unmasked stray light or clouds, and no sign of iterative failure 
of atmospheric correction (e.g. Figs. 2,4,5). 

The primary purpose of this study was the investigation of spatial 
statistics leading to a better understanding of instrument requirements 
for future satellite missions. Analytically, the results were more depen- 
dent upon interpixel variability in ocean color rather than absolute 
magnitudes, and no radiometric validation of retrievals was attempted. 
However, validation could precede as more field measurements of 
ocean color and biogeochemical properties become available in these 
extreme environments. Validation could then be extended to retrievals 
of biogeochemical properties such as TSM, provided field data were col- 
lected to regionally characterize absorption effects (e.g. from CDOM or 
phytoplankton pigments) on reflectance at 645 nm. With better charac- 
terization of on-orbit instrument SNR and a method for aerosol removal, 
it may be possible in the future to further refine the GSD approach 
using data from the Hyperspectral Imager for the Coastal Ocean (HICO) 
instrument aboard the International Space Station, which collects 


episodic data of the regions studied here at a spatial resolution of about 
100 m. The result would be to lower the optimal GSD lower limit to 
150 m, which, while not relevant in most of the nearshore, shelf, or off- 
shore waters studied here, may be important to inland waters and turbid 
river outflows. Similarly, given sensor-specific SNR as a function of L t , this 
analysis could be expanded to other visible wavebands using archived 
MER1S 300 m data, although the lower limit would increase somewhat 
from 375 m to 450 m. In addition to better understanding spatial variabil- 
ity in ocean color, the methods described here may be extrapolated to de- 
velop pixel-specific maps of ocean color uncertainty for remote sensing 
reflectance products derived from each of these instruments, provided 
the requisite sensitivity characterizations were available. 
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Fig. 9. Semivariograms and related parameters from HKM MODIS imagery at the two smaller river plume regions, the Chesapeake Bay/Susquehanna and Mississippi (see Figs. 7 & 8 
for descriptive details). 
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Fig. 10. Semivariograms and related parameters from HKM MOD1S imagery at the off- 
shore region in the Sargasso Sea (see Figs. 7 & 8for descriptive details). Note the change 
in sill scale between Figs. 7 and 8, and 9. 
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